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Abstract 

We perform a numerical study of the phase transitions in three-dimensional 
Z{N) lattice gauge theories at finite temperature for A^ > 4. Using the dual formu- 
lation of the models and a cluster algorithm we locate the position of the critical 
points and study the critical behavior across both phase transitions in details. In 
particular, we determine various critical indices, compute the average action and 
the specific heat. Our results are consistent with the two transitions being of infinite 
order. Furthermore, they belong to the universality class of two-dimensional Z{N) 
vector spin models. 



e-mail addresses: 

^oleg@bitp. kiev.ua, *vchelnokov@i.ua, ''"''cortese@unizar.es, ■'■gravina@ucy.ac.cy, 
'papa@cs.infn.it, **i_van_go@inbox.ru 



1 



1 Introduction 



The deconfinement phase transition in finite-temperature lattice gauge theories (LGTs) 
has been one of the main subjects of investigation for the last three decades. By now it is 
well studied and understood for a number of pure gauge models in dimensions D = 3,4. 
In particular, the phase structure of a finite-temperature three-dimensional {3D) pure 
SU{N) LGT with the standard Wilson action is thoroughly investigated both for N = 2,3 
and for the large- limit (see, e.g., [1] and references therein). The transition is second 
order for = 2, 3 and first order for > 4. In the case of the SU{4) gauge group, 
most works agree that the transition is weakly first order. The deconfining transition in 
SU{N = 2, 3) LGTs belongs to the universality class of 2D Z{N = 2, 3) Potts models. 
All these phase transitions are characterized by the spontaneous symmetry breaking of a 
Z{N) global symmetry of the lattice action in the high-temperature deconfining phase. 

Another interesting set of models is represented by abelian Z{N) LGTs. Z{N) is the 
center subgroup of SU{N), hence Z{N) LGT can provide useful insights into the universal 
properties of SU{N) models. Moreover, Z{N) LGTs are interesting on their own right 
and might possess an even richer phase structure as will be revealed below. The most 
general action for the Z{N) LGT can be written as 

'S'gaugo = X] X] '^l^kCOs{^^{Sn{x) + Sm{x + Cn) - Sn{x + Cm) - Sm{x))\ . (1) 
X n<m k=l ' 

Gauge fields are defined on links of the lattice and take on values s„(a;) = 0, 1, ■ ■ ■ , A^ — 1. 
Z{N) gauge models, similarly to their spin cousins, can generally be divided into two 
classes - the standard Potts models and the vector models. The standard gauge Potts 
model corresponds to the choice when all Pk are equal. Then, the sum over in ([T]) reduces 
to a delta-function on the Z{N) group. The conventional vector model corresponds to 
f3k = for all k > 1. For A^ = 2, 3 the Potts and vector models are equivalent. 

The study of the phase structure of Z{N) LGTs at zero temperature has a long history. 
While the phase structure of the general model defined by ([1]) remains unknown, it is well 
established that the Potts models and vector models with only /3i non- vanishing have one 
phase transition from a confining phase to a phase with vanishing string tension [21 |2l S] . 
Via duality, Z{N) gauge models can be exactly related to 3D Z{N) spin models. In 
particular, a Potts gauge theory is mapped to a Potts spin model, and such a relation 
allows to establish the order of the phase transition. Hence, Potts LGTs with N = 2 have 
second order phase transition, while for A^ > 3 one finds a first order phase transition. 
Vector models have been studied numerically in [3] up to A^ = 20. It was confirmed 
that the zero-temperature models possess a single phase transition which disappears in 
the limit A^ — )■ oo. Thus, the U{1) LGT has a single confined phase in agreement with 
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theoretical results [B]. A scaling formula for the critical coupling with N had also been 
proposed in [5]. We are not aware, however, of any detailed study of the critical behavior 
of the vector models with > 4 in the vicinity of this single phase transition. 

The deconfinement phase transition at finite temperature is well understood and stud- 
ied for = 2, 3. An especially detailed study was performed on the gauge Ising model, 
A^ = 2, in [7]. These models belong to the universality class of 2D Z{N) spin models 
and exhibit a second order phase transition in agreement with the Svetitsky-Yaffe con- 
jecture [8j. One should expect on general grounds that the gauge Potts models possess a 
first order phase transition for all A^ > 4, similarly to 2D Potts models. The 2'(4) vector 
model has been simulated, e.g., in [9]. It also belongs to the universality class of the 2D 
Z{4) spin model and exhibits a second order transition. 

Much less is known about the finite-temperature deconfinement transition for the 
vector Z{N) LGTs when A^ > 4. It is the ultimate goal of the present work to deepen 
our understanding of the phase structure of these models. The Svetitsky-Yaffe conjecture 
is known to connect critical properties of 3D Z{N) LGTs at finite temperature with 
the corresponding properties of 2D spin models, if they share the same global symmetry 
of the action. It is widely expected, and in many cases proved by either analytical or 
numerical methods, that some 2D Z{N > 4) spin models (like the vector Potts model) 
have two phase transitions of infinite order, known as the Berezinskii-Kosterlitz-Thouless 
(BKT) phase transitions. According to the conjecture, the phase transitions in some 
3D Z{N > 4) gauge models at finite temperature could exhibit two phase transitions as 
well. Moreover, if the correlation length diverges when approaching the critical point, 
these transitions should be of the BKT type and belong to the universality class of the 
corresponding 2D Z[N) spin models. 

The BKT phase transition has been best studied in the 2D XY model [TOj [TH [12] . 
Certain analytical [H [131 [H] as well as numerical results [151 [IS] unambiguously indicate 
that the deconfining phase transition in the 3D U{1) LOT is also of infinite order and 
might belong to the universality class of 2D XY model 0. 

In recent papers [171 H] we have initiated exploring the phase structure of the vector 
Z{N) LGT for A^ > 4. More precisely, we have considered an anisotropic lattice in the 
limit where the spatial coupling vanishes. In this limit the spatial gauge fields can be 
exactly integrated out and one gets a 2D generalized Z{N) model. The Polyakov loops 
play the role of Z{N) spins in this model. For the Villain version of the model obtained 
we have been able to present renormalization group arguments indicating the existence of 
two BKT-like phase transitions. This scenario was confirmed with the help of large-scale 
Monte Carlo simulations of the effective model. We have also computed some critical 

^It should be noted, however, that the numerical results of |16j point to a critical index ?/ larger than 
its XY value by almost a factor of 2 for Nt = 8. Therefore, the question of the universality remains open 
for this model. 
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indices which appear to agree with the corresponding indices of 2D Z{N) spin models, 
thus giving further support to the Svetitsky-Yaffe conjecture. 

In this paper we extend our analysis to the full isotropic 3D Z{N) LGT at finite 
temperature. It is well known that the full phase structure of a finite-temperature LGT 
is correctly reproduced in the limit where spatial plaquettes are neglected. They have 
probably small influence on the dynamics of the Polyakov loop interaction. We therefore 
expect that the scenario advocated by us in [T7] remains qualitatively correct for the full 
theory. In particular, full gauge models with > 4 should possess two phase transitions 
of the BKT type and we expect the values of critical indices to coincide with the indices 
of the 2D vector spin models. 

The 2D Z{N) spin model in the Villain formulation has been studied analytically in 
Refs. [H dSl [201 [221 U3\ ■ It was shown that the model has at least two phase transitions 
when N > 5. The intermediate phase is a massless phase with power-like decay of the 
correlation function. It turns out that //(/^i^^) = 1/4 at the transition point from the strong 
coupling high-temperature phase to the massless phase, i.e. the behavior is similar to that 
of the XY model. At the transition point /3c from the massless phase to the ordered 
low-temperature phase one has r]{l3P) = A/N\ A rigorous proof that the BKT phase 
transition does take place, and so that the massless phase exists, has been constructed 
in Ref. [23] for both Villain and standard formulations of the vector model. Universality 
properties of vector models were studied via Monte Carlo simulations in Ref. [25] for 

= 6, 8, 12 and in Refs. [23 [23 [211 [29] for A^ = 5, 7, 17. Results for the critical indices t] 
and u agree with analytical predictions obtained for the Villain formulation of the model. 

A similar phase structure is expected to hold for the finite-temperature 3D Z{N > 4) 
LGT. It can be described in terms of the Polyakov loop correlation functions as follows. 
The low-temperature phase is a confining phase. Here, the correlation decays with an 
area law, thus implying a non-vanishing string tension and a linear potential between 
static Z{N) charges. With the temperature increasing, the system undergoes a phase 
transition to a massless phase. This phase is characterized by the enhancement of the 
Z{N) global symmetry to a U{1) symmetry and is very close in nature to the high- 
temperature phase of U{1) gauge model. In particular, the dominating contribution to 
the correlation of the Polyakov loops comes from massless excitations (dual spin- waves). 
This is nevertheless a confinement phase, as the correlation decays with a power law and 
though the string tension vanishes, the potential between Z{N) charges is logarithmic. 
Increasing the temperature further on leads to a spontaneous breaking of the Z{N) global 
symmetry at some critical point. One enters a deconfining phase above this critical point. 

Here we intend to: 

• check the scenario of two BKT phase transitions described above; 

• compute some critical indices at both transitions and verify the universality class of 
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the model. 



The fact that the BKT transition has infinite order makes it hard to study its properties 
using analytical methods. In most of the cases studied so far one uses a renormalization 
group (RG) technique as in Ref. [19]. Unfortunately, there are no direct ways to generalize 
the transformations of Ref. |T9], leading to RG equations, to 3-D Z{N) LGTs, except 
for the limiting case — > oo. To study the phase structure of these models we need 
numerical simulations. Here, however, another problem appears related to the very slow, 
logarithmic convergence to the thermodynamic limit in the vicinity of the BKT transition. 
It is thus necessary to use both large-scale simulations and combine them with the finite- 
size scaling methods. Our principal strategy consists in passing to a dual formulation of 
the 3D Z{N) vector LGT, which is known to be a generalized 3D Z{N) vector model. 
This allows us to use a cluster algorithm in our simulations. The standard procedure in 
studying the phase structure is to use Binder cumulants to locate the position of critical 
points. Then, critical indices can be determined from various susceptibilities. In the 
case of a deconfinement phase transition both Binder cumulants and susceptibilities are 
usually constructed from Polyakov loops. However, it is a non-trivial problem to write 
down the expression for a single Polyakov loop in a dual formulation (though the dual 
form can be easily found for invariant quantities, like the correlation of Polyakov loops). 
We have therefore decided to study the critical behavior making use of Binder cumulants 
and susceptibilities constructed from the dual Z(N) spins d. This procedure exhibits 
an interesting phenomenon, namely the critical behavior of dual spins is reversed with 
respect to the critical behavior of Polyakov loops: the spontaneously-broken ordered phase 
is mapped to the symmetric phase and vice versa. Moreover, the critical indices t] are also 
interchanged as will be explained later on. The index u which governs the exponential 
divergence of the correlation length is expected to be the same at both transitions and 
takes on the value u = 1/2. 

The lowest number of N where two BKT phase transitions are expected is = 5. 
In the present paper we study the phase transitions in models with = 5, 13 in great 
detail. We have studied these values of in the strong coupling regime [HI HE], therefore 
it is natural to continue working with these models. Our computations are performed on 
lattices with temporal extent Nt = 2,4 and with spatial size in the range L G [32 — 1024]. 

This paper is organized as follows. In Section 2 we formulate our model and establish 
the exact relation with a generalized 3D Z{N) spin model. Here, we also explain why 
critical indices, determined from the dual spin correlation functions, become 77 (/^i^-*) = 
A/N'^ and ?7(/3i^^) = 1/4, i.e. the values are interchanged with the respect to what expected 

^The use of the Polyakov loop correlations to directly extract critical indices requires, in case of an 
essential singularity occurring at the BKT transition, prohibitively huge lattices, which are not accessible 
with the present numerical facilities. 
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from the correlation of original degrees of freedom. In Section 3 we present the setup 
of Monte Carlo simulations, define the observables used in this work and present the 
numerical results of simulations. In particular, we locate the position of critical points 
and compute various critical indices at these points. As a further cross-check of the nature 
of the phase transitions, we also present results of a few simulations for the Z(5) gauge 
model and compute the average action and the specific heat in the vicinity of critical 
points. Our conclusions and perspectives are given in Section 5. 

2 Relation of the 3D Z{N) LGT to a generalized 3D 
Z[N) spin model 

2.1 Partition and correlation functions 

We work on a ?>D lattice A — L"^ x Nt with spatial extension L and temporal extension 

A^^; X = {xq,Xi,X2), where Xq G [0, iV^ — 1] and X1.X2 G [0, L — 1] denote the sites of the 
lattice and e„, n = 0, 1, 2, denotes a unit vector in the n-th direction. Periodic boundary 
conditions (BC) on gauge fields are imposed in all directions. The notations pt (ps) stand 
for the temporal (spatial) plaquettes. It (h) for the temporal (spatial) links. 
We introduce conventional plaquette angles s{p) as 

sip) = Sn{x) + Sm{x + en) - Sn{x + Cm) - Smix) . (2) 

The 3D Z{N) gauge theory on an anisotropic lattice can generally be defined as 

Z(A;A,/3,;iV) = J] f ^ E ) n^(^(^'^)) 11 ^(^(^'*)) " (3) 

leA \ s{i)=o/ Ps pt 



The most general Z(A^)-invariant Boltzmann weight with A^ — 1 different couplings is 



Q{s) = exp 



2nk 



cos s 

N 
fe=i 



(4) 



The Wilson action corresponds to the choice /3p(l) = l3p{k) = 0,k = 2, A" — 1. The 
C/(l) gauge model is defined as the hmit A" — > oo of the above expressions. 

In what follows we work with the conventional Wilson action on an isotropic lattice, 
i.e. Ps = Pt = P- To study the phase structure of 3D Z{N) LGTs one can map the gauge 
model to a generalized 3D spin Z{N) model with the action 

^ = E E E iv" ^'^"^^ " '^"^ + '"^0 ■ 

X n=l k=l 
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The effective couphng constants /3k can be computed exactly as follows. The first step is 
to construct a dual form for the partition function Details can be found, e.g., in 
One gets on the dual lattice the following expression for the partition function |fl 




Yl Qdis{x) - s{x + e„)) 



(6) 



where the dual Boltzmann weight Qd{s) becomes 



N-1 



Qd{s) 



^ Int+sW) = J^exp 



p=0 



(3 cos 



2ttp 
1^ 



( 2TTpS 



(7) 



Here, /^(x) is the modified Bessel function. Exponentiating and re-expanding the dual 
weight in a new Fourier series one finds Pk as 



(3k 



N-1 

-y\n 



p=0 



Qd{0) 



( 2TTpk 



As example, we give below the expressions for the effective couplings Pk for A^ 
obtains from ([7]) and ([8]) 



(8) 



5. One 



Pi 

P2 



Pa 

Ps 



^ ln[l - U/2] - ln[l - t_/2] 



10 '''' " ^ 10 
^ + ln[l - t+/2] - ^—y^ ln[l - t_/2] 



10 



10 



(9) 



where 



5±V5+ (5^^5)6 2 

y/^ f. 



f. 



2 + 2e^^ + ei(5+^)^ 

As is seen from Eq. ([7]), the dual model is ferromagnetic. However, P2 is small and 
negative for P > 1.077. In the whole region \Pi\ ^ \P2\; moreover, in the critical region 
|/^2|/|/3i| ~ 10~^. Similar properties hold for all A^. Thus, one expects that the 3D vector 
spin model with only Pi non-vanishing gives a reasonable approximation to the gauge 
model (in our simulations we use all Pk). Next important fact, evident from Eq. ([9]), 



^In writing this expression we have neglected a certain global summation which appears due to the 
fact that the product of the plaquette variables around the closed 2d surface winding through the lattice 
in periodic directions is unity (so-called global Bianchi identity). Such global identities can be safely 
omitted since they do not influence the quantities of our interest in thermodynamic limit. Note, this is 
not the case for quantities like the twist free energy. 
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is that the weak and the strong couphng regimes are interchanged: when /3 — )■ oo both 
effective couplings (3^ ^ and, therefore, the ordered symmetry-broken phase is mapped 
to a symmetric phase with vanishing magnetization of dual spins. The symmetric phase 
at small /3 becomes an ordered phase where the dual magnetization is non-zero. 

Let W{x) = n^*=o 6^P(^i ^i.^)) be the Polyakov loop in the representation j. The 
correlation function of Polyakov loops can be expressed in the dual form as 



where 5*^ is the dual surface enclosed between two Polyakov loops, i. e. it consists of links 
perpendicular to plaquettes of the original lattice and is closed in the temporal direction. 
The correlation function of the dual spins is defined in standard way as 



In the original (gauge) formulation this correlation function corresponds to a disorder 
operator which is the 't Hooft line in the 3D theory. It consists of a string of plaquettes 
connecting the cubes dual to the points and R and measures the free energy of a Z(N) 
monopole-anti-monopole pair. 

2.2 Behavior of the dual correlation function 

As was explained in the Introduction, our aim is to study the critical behavior making 
use of quantities, like the Binder cumulants, constructed from the dual Z{N) spins. On 
very general grounds one expects that in the confined phase, where the potential between 
electric Z{N) charges grows linearly with the distance, the correlation function (fT2!) shows 
ordered behavior, so that the magnetic charges are free. In the phase where the electric 
charges are deconfined, the correlation function 0121) exhibits exponential decay and the 
free energy of a monopole-anti-monopole pair grows with distance. In the massless phase, 
if such phase exists, one anticipates that both kind of charges are confined by a weakly- 
growing logarithmic potential, hence both correlations decrease with a power law. It is not 




Here we have introduced the sources h{l) as 




(11) 




(12) 
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obvious, however, and to the best of our knowledge not proved so far, that the correlation 
of the Polyakov loops and the dual correlation show similar critical behavior. In fact, as 
we shall argue below, their critical behavior is exactly opposite. Let us suppose, there 
exists a massless phase and /3c^^ is the phase transition point from the electric confinement 

(1') 

to the massless phase, while /3c is the phase transition point from the massless phase 
to the deconfined phase. One then has for the Polyakov loop correlation the following 
asymptotic behavior at the critical points (for j = 1) 



P,-(/?;/3«;iV) 



P,(i?;/3f ;iV) 



' - 'j\-^->i^c 1 — y " ^^(2) 
In the case of the dual correlations, the rj indices are interchanged 

1 „ „.o^ 1 



r,(i?;/3«;iV) 



(2) ■> 



r,(/?;/3f ;iV) 



(1) 



(13) 



(14) 



If the ?)D Z{N > 4) vector gauge model belongs to the universality class of the corre- 
sponding 2D Z{N) spin model, we should find 



V 



(1) 



1/4 , f] 



(2) 



(15) 



For fixed and Nt increasing, Pc^^ — /3i^\ In the limit Nf ^ oo one ends up with a 
single phase transition from the confining to the deconfining phase. When A^^^ is fixed and 

increases, /3c^^ diverges roughly as A^^, while /3c^^ approaches the critical point of the 
finite-temperature U{1) LGT exponentially fast. 

The direct proof of all these properties would include the construction of RG equations 
describing the behavior of the system in the vicinity of phase transitions, something we 
could not accomplish so far. In the absence of such a proof, we give qualitative arguments 
why the behavior f ll3p - (IT^ is rather natural and might be anticipated. First of all, the 
interchange of the critical indices can be easily seen in the 2D Z{N) spin models. Indeed, 
Eq. (|T0|) remains formally correct for the two-point correlation function of these models if 
one replaces the surface Sd in (|TT|) with a dual path connecting the points and R. The 
two-point correlation of dual spins can be written in precisely the same form of Eq. f lTU]) in 
the original formulation, with Qd{s) substituted by Q{s). Let us consider now the Villain 
formulation. For Q{s) it means 



Q{s) ^ Q^(s) 



exp 



(16) 



while for Qd{s) it amounts to replacing the Bessel function in ([7]) with its asymptotics 



E 



exp 



m=— oo 



■ —(s + Nmf 



(17) 
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In this formulation the 2D Z{N) model is self-dual. From this fact, from Eq. f llUp and its 
analog for the dual correlation function it follows that (we keep the notation Pj for the 
correlation function of spins) 



P,iR;(3;N) 



47r2/3 



;iV 



It is now straightforward to repeat the calculations of Ref. [19] for the dual correlation 
function and to get 



Pj{R;(3]N) X exp 



J' 



2n(3, 



InR 



'eff 



R; 



47r2/3 



exp 



InR 



■ (19) 



/3eff can be expanded in powers of the self-energies of the topological defects of the model. 
The leading contribution P^s ~ P comes from the spin-wave configurations. At the critical 
points one has 

= 2/7r , /3efr(/3f ) = iVV(8vr) . (20) 

It then follows from Eqs. (1T^ - (I2U]) that the critical indices are indeed interchanged. 2D 
Z{N) vector models in the standard formulation are not self-dual, therefore Eq. (llSp 
does not generally hold. Nevertheless, this pictures remains qualitatively correct for 
all formulations of 2D models which belong to the same universality class. Moreover, 
recalling that \(3i\ ^ \(3k\ for all k ^ 1( mod A^), one finds the approximate equality 
Pj{R; (3; N) ^ f^ii ^) for the standard formulation, where /3k are the dual coupling 

constants (IE])- 

Let us turn now to the 3D Z{N) LGT and use again the Villain formulations, Eqs. ( |T6|) 
and (fT7|) (in the case of the gauge model one has to take the plaquette angle s{p) defined 
in Eq. ([2])). In the spin- wave approximation, which should be valid in the massless phase, 
we obtain 



Pj{R;(3;N) x exp 
for the Polyakov loops correlation function and 



T,{R;(3;N) 



exp 



'1^ 



f G{R) 



(21) 



(22) 



for the dual correlations. D{R) in Eq. ( l2T]) is the two-dimensional Green's function. 
G{R) in Eq. (1221) is the three-dimensional Green's function. At finite temperature in the 
massless phase the dominant contribution to G{R) arises from the temporal zero mode, 
other modes being massive and exponentially suppressed. We thus have 



GiR) ^ i- DiR) + 0{e 



-mR\ 



D(R) ^ —InR. 
2n 



10 



Combining this with Eq. (1^ and Eq. ( 12^ leads to 



Pj{R;(3;N) x exp 



\nR 



TjiR;(3;N) 



exp 



2n(3 



J 



(23) 



These expressions must be supplemented by computing the corrections from the topolog- 
ical defects of the model, like Z{N) monopoles. We expect that this results, similarly to 
2D models, in replacing /3 — )■ in the above equations. It is natural to suppose that 
the analog of Eq. (120|) takes the form 



N, 



TT 



(2)^ 



iV2 



57r 



(24) 



It corresponds to approximately linear scaling of the critical points with Nf. If so, the 
critical indices are interchanged as it becomes evident after the substitution of the last 
equations in fl23|l . As a matter of fact, Eq. flMl) is our conjecture which remains to be 
proved. 



3 Numerical results 

3.1 Setup of the Monte Carlo simulation 

To study the phase transitions we used the cluster algorithm described in [27]. We sim- 
ulate the dual model defined by Eq. on an Nt x L x L lattice with periodic BC. 
Simulations were carried out for Nt = 2, 4. To check directly that the critical indices are 
interchanged, as explained in the previous Section, we have also simulated the Z{5) gauge 
model for = 2, 4 using the heat-bath algorithm. As original action of the gauge model 
we used the conventional Wilson action. For each Monte Carlo run the typical number 
of generated configurations was 10^ the first 10^ of them being discarded to ensure ther- 
malization. Measurements were taken after 10 updatings and error bars were estimated 
by the jackknife method combined with binning. 
We considered the following observables: 

• complex magnetization = \Mi\e'^'^ , 

Ml = Y.^^p(^s{x)) ; (25) 



population Sl 
where rii is number of s{x) equal to i; 



N f maXi=o,7V-i ni l\ 
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Figure 1: Scatter plot of the complex magnetization Ml at /9=1.84, 1.98 and 2.06 in Z(5) 
on a 512^ x 4 lattice. 



real part of the rotated magnetization Mr = \Ml \ cos{NiIj) and normalized rotated 
magnetization = cos{Nip); 



susceptibilities of M^, Sl and Mr. x^l^\ X^l \ X^l^^^ 



(27) 



Binder cumulants f/l^^-* and B\ , 



U 



(M) 



B 



(Mr) 



1 - 



3 

\Mr 



M, 



R 



{Mn)\y 



{2i 



The variable s{x) appearing in the definitions fl25p and 0261) represents the dual spin in 
the case of simulations of the dual model. In the case of the gauge model, s{x) is the 
Polyakov loop, the sum in Eq. ( 125|) becomes two-dimensional and Nt must be omitted 
from Eqs. (l26|l and (1271) . Also, when simulating the gauge model we have computed the 
average action and the specific heat in the vicinity of the critical points. 



3.2 Determination of the critical couplings 

A clear indication of the three-phase structure emerges from the inspection of the scatter 
plot of the complex magnetization M^ at different values of /3: as we move from low to 
high /3, we observe the transition from an ordered phase [N isolated spots) through an 
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Figure 2: Binder cumulant U^^^^ as function of (3 (left) and of (/3 — /3c)(lnL)^/'^ (right) in 
Z{5) Nt = 2 model. 



intermediate phase (ring distribution) up to the disordered phase (uniform distribution 
around zero). Fig. [1] shows this three-phase structure for the case of Z{5) on a 512^ x 4 
lattice. 

The first and most important numerical task is to determine the value of the two critical 
couplings in the thermodynamic limit, pi^^ and I3c^\ that separate the three phases. To 
this aim we find the value of Pc which provides the best overlap of universal observables, 
plotted for different values of L against {/3 — /3c^-')(ln L)^/*^, with u fixed at 1/2. As these 
universal observables we used: 

• Binder cumulant B^^^'^^ and the order parameter for the first phase transition; 

• Binder cumulant U^^^^ for the second phase transition. 

To localize regions in which the overlap must occur, the approximate maxima of the 
susceptibilities x^l^ and Xl^"* were used. 

In Figures [2lll] we give the plots of the universal observables, drawn against /3 and 
against {P—Pc^^){\n L)^/" . We report in Table[T]the determinations of the critical couplings 
(3^c^ and pf^ in Z{N) with A^=5 and 13 for Nt=2 and 4. 

3.3 Determination of critical indices at the two transitions 

Once critical couplings have been estimated, we are able to extract some critical indices 
and check the hyperscaling relation. 

Since we are using the observables in the dual model the transitions change places: 
the first transition is governed by the behavior of Mr, the second one by the behavior of 
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Table 1: Values of and obtained for Nt = 2 and 4 in Z{N) with = 5, 13. 





Nt 






5 


2 


1.617(2) 


1.694(2) 


5 


4 


1.943(2) 


1.990(2) 


13 


2 


1.795(4) 


9.699(6) 


13 


4 


2.74(5) 


11.966(7) 



Ml. This could already be seen in the previous Section, since the first critical point was 
obtained from the i?4*^^^ and the second one from the U^^^ curve collapse. 

We start the discussion from the second transition. According to the standard finite- 
size scaling (FSS) theory, the equilibrium magnetization \Ml\ at criticality should obey 
the relation \Ml\ ~ L~^^'^, if the spatial extension L of the lattice is large enough q 
Therefore, we fit data of \Ml\ at /3c , on all lattices with size L not smaller than a given 
-^^min, with the scaling law 

\Ml\= AL-^/" . (29) 

The FSS behavior of the susceptibility Xl^"* is given by Xl^'* ~ L"'^'^, where 'y/jy = 2—7] 
and 7] is the magnetic critical index. Therefore we fit data of X/f^^ \ on all lattices 

with size L not smaller than a given Lmin, according to the scaling law 

X^L ^ = Ar/^ . (30) 

As the value of the critical coupling /3c we use the central value determined in the 
previous Section. 

The results of the fits are summarized in Tabled Each row corresponds to the fit using 
data from L = Lmin up to L = 1024. The reference value for the index rj at this transition 
is 1/4, whereas the the hyperscaling relation to be fulfilled is7/z/ + 2/3/z/ = (i = 2. 

We see that in most cases the values of r] and d are close to those predicted by 
universality. The discrepancy from the exact values t] = 0.25 and d = 2 may be caused 
by the asymptotically vanishing parts of the scaling behavior of the observables \Ml\ and 
X^L^\ that we are not taking into account, but may be significant for smaller lattice sizes. 

The procedure for the determination of the critical indices at the first transition is 
similar to the one for the second transition, with the difference that the fit with the 



^Thc symbol /3 here denotes a critical index and not, obviously, the coupling of the theory. In spite 
of this inconvenient notation, we are confident that no confusion will arise, since it will be always clear 
from the context which /3 is to be referred to. 
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Table 2: Critical indices fi/v and 7/z/ for the second transition in Z{N) models, deter- 
mined by the fits given in Eqs. fl29p and fl3UI) on the complex magnetization Ml and its 
susceptibility Xl^'' ^it /^F"* for different choices of the minimum lattice size Lmin- The 
of the two fits, given in the columns four and six, is the reduced one. 



model 


Lm'm 


/3/z/ 




7/z/ 




d = 2(5/v + -i/v 


= 2 — 7/z/ 




128 


0.1226(4) 


0.92 


1.76(1) 


1.75 


2.00(1) 


0.24(1) 


Z{5) 


192 


0.1226(6) 


1.15 


1.75(2) 


2.14 


2.00(2) 


0.25(2) 


Nt = 2 


256 


0.1226(9) 


1.53 


1.77(2) 


1.51 


2.02(2) 


0.23(2) 


= 1.694 


384 


0.121(1) 


0.93 


1.74(2) 


0.78 


1.98(3) 


0.26(2) 




512 


0.1230(2) 


0.011 


1.74(5) 


1.46 


1.99(5) 


0.26(5) 




32 


0.1078(2) 


5.17 


1.69(1) 


51.5 


1.91(1) 


0.31(1) 




64 


0.1075(1) 


1.56 


1.71(1) 


19.6 


1.93(1) 


0.29(1) 


Z{5) 


128 


0.1074(1) 


1.39 


1.734(8) 


5.50 


1.949(8) 


0.266(7) 




192 


0.1075(2) 


1.67 


1.744(7) 


2.65 


1.959(7) 


0.256(7) 


= 1.990 


256 


0.1077(2) 


1.27 


1.752(6) 


1.40 


1.968(7) 


0.248(6) 




384 


0.1076(4) 


1.67 


1.760(8) 


1.10 


1.975(9) 


0.240(8) 




512 


0.1081(5) 


1.15 


1.77(1) 


1.57 


1.98(1) 


0.23(1) 




128 


0.1225(4) 


1.03 


1.749(8) 


0.428 


1.994(9) 


0.251(8) 


Z(13) 


192 


0.1225(6) 


1.29 


1.758(7) 


0.232 


2.003(8) 


0.242(7) 


Nt = 2 


256 


0.1229(8) 


1.49 


1.749(6) 


0.126 


1.995(8) 


0.251(6) 


pP = 9.699 


384 


0.123(1) 


2.20 


1.741(9) 


0.103 


1.99(1) 


0.259(9) 




512 


0.1203(2) 


0.0187 


1.73(1) 


0.0794 


1.97(1) 


0.27(1) 




32 


0.1266(5) 


16.62 


1.70(1) 


19.96 


1.95(1) 


0.30(1) 


Z(13) 


384 


0.1275(4) 


5.14 


1.72(1) 


7.70 


1.98(1) 


0.28(1) 


Nt = A 


128 


0.1282(2) 


0.54 


1.747(8) 


1.58 


2.004(8) 


0.253(8) 


13^ = 11.966 


192 


0.1283(3) 


0.76 


1.76(1) 


1.54 


2.01(1) 


0.24(1) 




256 


0.1282(6) 


1.45 


1.74(2) 


1.14 


2.00(2) 


0.26(2) 
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Table 3: Critical indices P/v and 7/z/ for the first transition in Z{N) models, deter- 
mined by the fits given in Eqs. (^U^ and fl5Ul) on the rotated magnetization Mji and its 
susceptibility Xzf^^'' (^c^^ ^'^^ different choices of the minimum lattice size Lmin. 



model 




/3/u 




7/z/ 


X-y/u 


d = 2/3/z/ + 7/z/ 


?7 = 2 — 7/z/ 




128 


0.097(6) 


0.101 


1.847(5) 


0.561 


2.04(2) 


0.153(5) 


Z(5) 


192 


0.103(8) 


0.093 


1.841(6) 


0.447 


2.05(2) 


0.159(7) 


Nt = 2 


256 


0.10(1) 


0.122 


1.850(2) 


0.038 


2.06(3) 


0.150(2) 


f3^c> = 1.617 


384 


0.09(2) 


0.117 


1.851(4) 


0.056 


2.03(4) 


0.149(4) 




512 


0.10(3) 


0.198 


1.848(7) 


0.091 


2.05(8) 


0.152(7) 




32 


0.123(6) 


1.08 


1.8403(8) 


0.72 


2.09(1) 


0.1596(8) 




64 


0.118(9) 


1.13 


1.841(lj 


0.58 


2.08(2) 


0.159(1) 


Z{5) 


128 


0.11(1) 


1.25 


1.841(1) 


0.70 


2.07(3) 


0.159(1) 




192 


0.11(2) 


1.56 


1.842(2) 


0.86 


2.07(4) 


0.158(2) 


= 1.943 


256 


0.11(3) 


2.07 


1.842(3) 


1.14 


2.06(6) 


0.158(3) 




384 


0.07(4) 


1.68 


1.836(2) 


0.22 


1.98(8) 


0.164(2) 




512 


0.06(8) 


3.26 


1.837(4) 


0.42 


2.0(2) 


0.163(4) 




128 


0.07(5) 


1.28 


1.968(9) 


0.97 


2.1(1) 


0.032(9) 


Z(13) 


192 


0.02(5) 


1.16 


1.97(1) 


1.11 


2.0(1) 


0.03(1) 


Nt = 2 


256 


0.04(8) 


1.48 


1.97(2) 


1.43 


2.0(2) 


0.03(2) 


13^^^ = 1.795 


384 


0.06(14) 


2.19 


1.98(3) 


1.93 


2.1(3) 


0.02(3) 




512 


0.2(17) 


1.08 


1.94(5) 


1.93 


2(3) 


0.06(5) 


Z(13), iVi=4 


128 


-0.3(1) 


0.59 


1.977(2) 


0.21 


1.3(3) 


0.023(2) 


= 2.74 


192 


-0.5(1) 


0.20 


1.980(4) 


0.24 


0.9(2) 


0.020(4) 



scaling laws Eqs. fl29|) and f l30|) is to be applied to data of the rotated magnetization, M/j, 
and of its susceptibility, x^l '"^ 1 respectively. As the value of the critical coupling /^c^-* we 
use the central value determined in the previous Section. 

The results of the fits are summarized in Table [3l The reference value for the index r] 
at this transition is A/N"^, i.e. rj = 0.16 for N = 5 and rj ^ 0.0237 for = 13, whereas 
the hyperscaling relation to be fulfilled is 7/z/ + 2/3/z/ = d = 2. 

A general comment is that, in many of the cases we investigated, both d and 1] at the 
two critical points slightly differ from the expected values, though these differences cancel 
to a large extent if we define rj as 2/3/z/. 

Also here we see a general agreement between the 77 and d values obtained and those 
predicted by universality. However, the expected value of /3/z/ in Eq. ( 129|) is very small. 
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Figure 5: Correlation between Xl L'^ ^ and the Binder cumulant i?4 in Z{13) with 
A^t = 4 for 1] = 0.25 (left) and for r] = 0.0237 (right) on lattices with different size. 

{2/N'^), so other, asymptotically vanishing, terms can have a great impact on its de- 
termination on finite-sized lattices. This is especially evident for 2^(13) with A^^^ = 4, 
where /3/i/ is negative indicating that the magnetization Mfj grows with lattice size. This 
makes problematic the determination of this index from simulations on lattices used in 
the present work and explains the discrepancies with the d = 2 value. 

There is an independent method to determine the critical exponent 77, which does not 
rely on the prior knowledge of the critical coupling, but is based on the construction of 
a suitable universal quantity [301 E?]. The idea is to plot Xl!'^^^L^~'^ versus and 
to look for the value of 77 which optimizes the overlap of curves from different volumes. 
This method is illustrated in Fig|5l for Z{13) model with Nt = 4. Another option is to 
plot MrL"/^ versus m^, which leads to overlapping curves for rj fixed at the value of the 
second phase transition, as illustrated in Fig |6] for 2^(13) model with Nt = 2. 

Concerning the value of the critical index u, the methods used in this work do not 
allow for the direct determination of its value. When locating critical points we have fixed 
at 1/2. This value appears to be well in agreement with all numerical data. 

3.4 Other checks of the nature of the phase transitions 

In this Section we describe briefly some results of the simulation of the original gauge 
model. We have simulated Z{5) LGT with A^t = 2,4 and spatial extent L E [64 — 512]. 
The typical number of measurements was 10^. In Tables H] and |5] we present results for 
Nf = 2 and Nt = 4, correspondingly. In general, errors are bigger and results for critical 
indices are not so precise as in dual model simulations. Nevertheless, we can conclude 
from the inspection of Tables H] and |5] that (i) the critical index rj is compatible with its 
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lattices with various values of L. 



Table 4: Critical indices /3/z/ and 7/1^ for the transitions in Z{b) LGT with Nt = 2. 







f3/u 




7/z/ 




d = 2/3/i/ + 7/i/ 


?7 = 2 — 7/z/ 


Pl^^ = 1.617 


192 


0.127(2) 


2.38 


1.78(5) 


3.15 


2.03(5) 


0.22(5) 


/3f ^ = 1.694 


192 


0.1(2) 


7.62 


1.82(6) 


7.51 


2.1(4) 


0.18(6) 



2D value and (ii) the values of the indices at two transitions are indeed interchanged as 
explained in Section 2. 

To produce further evidence in favor of the fact that the phase transitions investigated 
so far are both of infinite order, we have calculated the average action and the specific heat 
around the transitions in Z{5) LGT with Nt=2 and Nt = 4. In all cases the dependence 
of these quantities on /3 is continuous. It follows that first and second order transitions are 
ruled out. As an example, we present in Fig. [7] the results of simulations for the average 
action for the case of Z{5) with A^^^ = 2. 
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Table 5: Critical indices P/v and 7/z/ for the transitions in ^(5) LGT with Nt = 4. 







(3/1^ 




7/z/ 




= 2/3/z/ + 7/z/ 


= 2 — 7/z/ 


/3i^^ = 1.943 


128 


0.122(4) 


3.30 


1.71(7) 


2.55 


1.95(8) 


0.29(7) 


I3P = 1.990 


64 


0.4(2) 


4.92 


1.82(2) 


1.25 


2.6(5) 


0.18(2) 




4 Conclusions and Perspectives 



In this paper we have studied the 3D Z{N) LGT at the finite temperature aiming at 
shedding hght on the nature of phase transitions in these models for N > A. This 
study was based on the exact dual transformation of the gauge models to generahzed 3D 
Z{N) spin models. In Section 2 wc presented an overview of the exact relation between 
couphngs of these two models and described quahtatively the behavior of the Polyakov 
loop correlation function and correlation of dual spins (the disorder operator in the gauge 
formulation). Furthermore, we have advanced some arguments that the critical behavior 
of Polyakov loop correlations and dual correlations is reversed, in particular the values 
of the critical index rj at the two transitions are interchanged for the dual correlation 
function with respect to the Polyakov loop correlation. 

The numerical part of the work has been devoted to the localization of the critical 
couplings and to the computation of the critical indices. The main results can be shortly 
summarized as follows: 

• We have determined numerically the two critical couplings of Z{N = 5, 13) LGTs 
and given estimates of the critical indices r] at both transitions. For the first time we 
have a clear indication that for full 3D Z{N) vector LGT with N > 5 the scenario 
of three phases is realized: a disordered phase at small /3, a massless or BKT one at 

intermediate values of /3 and an ordered phase, occurring at larger and larger values 
of (3 as N increases. This matches perfectly with the N ^ oo limit, i.e. the 3D 
U{1) LGT, where the ordered phase is absent; 

• We have found that the values of the critical index r) at the two transitions are 
compatible with the theoretical expectations; 

• The index i/ also appears to be compatible with the value 1/2, in agreement with 
universality predictions. 

On the basis of this study we are led to conclude that finite-temperature 3D Z{N) 
LGTs for N > A undergo two phase transitions of the BKT type. Moreover, these models 
belong to the universality class of the 2D Z{N) vector models which also have two infinite 
order phase transitions and a massless phase. 

It should be emphasized that in this paper wc have not studied the continuum limit of 
the theory, but rather concentrated on the very possibility of having two BKT- like phase 
transitions. At the moment we are performing simulations of the model on symmetric 
lattices with the goal to compute the zero-temperature string tension and extend our 
present computations to Nt = 8. All these will help constructing the continuum limit. 
Also, we plan to extend our work to even N, to calculate the helicity modulus and to 
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establish scaling formulas for critical points with A^. The results of these studies will be 
reported elsewhere. 
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